#ifndef __UTILS__L2ERROR__
#define __UTILS__L2ERROR__

#include "../include/Grid.hpp"
#include "../include/Region.hpp"
#include <Eigen/Sparse>
#include <functional>
#include <cmath>
#include <array>

double L2error(Grid& g, Eigen::VectorXd x ,Region& r, std::function<double(double, double)> f)
{
    double rst = 0.0;
    int dof = r.basis_function.size();
    if(dof == 3)
    {
        for(auto iter=g.elements.begin(); iter<g.elements.end(); iter++)
        {
            double u0 = x[id((*iter)[0])];
            double u1 = x[id((*iter)[1])];
            double u2 = x[id((*iter)[2])];
            std::function<double(double, double)> f_g2 = [iter, u0, u1, u2, r, f](double x, double y)
            {
                std::array<double, 2> xi_eta = (*iter).g2l(x,y);
                double xi = xi_eta[0];
                double eta = xi_eta[1];
                double f_g = f(x,y) - u0*r.basis_function[0](xi,eta) - u1*r.basis_function[1](xi,eta) - u2*r.basis_function[2](xi,eta);
                return f_g * f_g;
            };
            rst += (*iter).integral(f_g2, r);
        }
    }
    if(dof == 6)
    {
        for(auto iter=g.elements.begin(); iter<g.elements.end(); iter++)
        {
            double u0 = x[id((*iter)[0])];
            double u1 = x[id((*iter)[1])];
            double u2 = x[id((*iter)[2])];
            double u3 = x[id((*iter)[3])];
            double u4 = x[id((*iter)[4])];
            double u5 = x[id((*iter)[5])];

            std::function<double(double, double)> f_g2 = [iter, u0, u1, u2, u3, u4, u5, r, f](double x, double y)
            {
                std::array<double, 2> xi_eta = (*iter).g2l(x,y);
                double xi = xi_eta[0];
                double eta = xi_eta[1];
                double f_g = f(x,y) - u0*r.basis_function[0](xi,eta) - \
                             u1*r.basis_function[1](xi,eta) - u2*r.basis_function[2](xi,eta) -\
                             u3*r.basis_function[3](xi,eta) - u4*r.basis_function[4](xi,eta) -\
                             u5*r.basis_function[5](xi,eta);
                return f_g * f_g;
            };
            rst += (*iter).integral(f_g2, r);
        }
    }
    return rst;
}

#endif
